RNA-Seq Data Analysis 2018

RAW READS MANIPULATION

Accendere il computer e aprire il terminale (ctrl + alt + T)

Spostarsi nella cartella Scaricati, create una nuova cartella chiamandola tc

cd Scaricati/
mkdir tc

Entrare nella cartella appena creata e valutarne il contenuto della cartella digitare il seguente comando

cd tc/
ls -lrth

La cartella dovrebbe essere vuota. Come prima cosa dobbiamo recuperare i file grezzi (FASTQ) del nostro esperimento. Per recuperare i file ci connettiamo al nostro server di laboratorio e con il protocollo di trasferimento ssh copiamo in locale (la cartella dove siamo in questo momento) i 2 file corrispondenti alle reads forward e alle reads reverse. Lavoreremo con un solo campione, ma come sapete l’esperimento consiste di 3 repliche biologiche per i non trattati e 3 repliche biologiche per quelli trattati. I comandi che qui eseguiremo possono essere ripetuti su più file contemporaneamente tramite appositi comandi come i cicli for

# prima scarichiamo il primo file (attenzione al punto in fondo al comando, indica current directory)
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/Sample_Untreated_1_1_subset.fastq.gz .
# inserire la password quando ce la chiede ()
# controlliamo che sia stato scaricato
ls -lrth
# ora scarichiamo il secondo file
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/Sample_Untreated_1_2_subset.fastq.gz .
# inserire nuovamente la password
# controlliamo
ls -lrth

Come potete osservare la dimensione di questi file è molto bassa rispetto alla dimensione tipica di un file che risulta da un esperimento di sequenziamento come lo abbiamo descritto. Questo perché i file originali sono stati copiati e ridotti in dimensioni con lo scopo di riuscire a manipolarli e lavorarci nei tempi della lezione, principalmente per il fatto che le capacità dei computer che stiamo usando potrebbero essere limitate e quindi non in grado di eleaborare questi file in tempi brevi

Iniziamo con esplorare questi file

# "apriamo" il file
less -S Sample_Untreated_1_1_subset.fastq.gz
# contiamo il numero di reads
zcat Sample_Untreated_1_1_subset.fastq.gz | echo $((`wc -l`/4))
zcat Sample_Untreated_1_2_subset.fastq.gz | echo $((`wc -l`/4))
# il numero di reads F e R coincide?

FASTQ file

QUALITY CONTROL

Per effettuare il controllo di qualità sulle nostre reads utilizzeremo il programma FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Questo software (scritto in Java) prende come input un numero a piacimento di reads e calcola varie statistiche che possono dare un’idea della qualità dei nostri dati

# creiamo una nuova cartella dove scaricheremo il programma e salveremo i risultati
mkdir fastqc
cd fastqc/
# scarichiamo il programma con il comando wget
wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.7.zip
# scompattiamo la cartella
unzip fastqc_v0.11.7.zip
# rimuoviamo cartella compressa
rm fastqc_v0.11.7.zip
# vediamo cosa c'è dentro la cartella senza spostarci dentro (percorso relativo alla nostra posizione)
ls -lrth FastQC/
# il file che ci interessa, cioè il programma che andremo a richiamare è "fastqc"
# dobbiamo prima renderlo eseguibile
chmod 755 FastQC/fastqc
# controlliamo che funzioni
./FastQC/fastqc -h

A questo punto possiamo lanciare il programma dando come input i due file con le reads F e R. Il comando è molto semplice, basta specificare i nomi dei file che vogliamo controllare, e la cartella dove vogliamo che vengano salvati i report coi risultati (la cartella deve già esistere quindi prima la creiamo)

# creiamo cartella per salvare i risultati
mkdir results
# ci spostiamo una cartella su, dove stanno i nostri file
cd ..
# eseguiamo il programma (dobbiamo specificare il percorso dell'eseguibile rispetto a dove siamo) e specifichiamo inoltre il percorso della cartella di output
./fastqc/FastQC/fastqc Sample_Untreated_1_1_subset.fastq.gz Sample_Untreated_1_2_subset.fastq.gz -o fastqc/results/

Per ogni file di input vengono creati 2 file di output, un file html e una cartella compressa. Andiamo ad esplorare il report in formato html aprendolo con Firefox

FastQC report

A questo punto commentiamo insieme i risultati del controllo qualità confrontandoci con qualità di riferimento sul sito di FastQC e con la spiegazione dei moduli di analisi (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3%20Analysis%20Modules/)

Phred quality

Dopo aver deciso insieme quali parametri utilizzare per correggere la qualità delle nostre reads ed aver identificato gli adattatori da rimuovere, andiamo ad effettuare la rimozione degli adattatori e il quality trimming

ADAPTERS REMOVAL AND QUALITY TRIMMING

Per effettuare queste due operazioni useremo un unico programma che si chiama Trimmomatic (http://www.usadellab.org/cms/?page=trimmomatic). Il suo utilizzo è molto simile a quello di FastQC (anche Trimmomatic è scritto in Java) e necessita come input le sequenze grezze. Per ogni file di input (nel nostro caso 2 file) verranno creati 2 file di output, uno indicato come paired, l’altro come unpaired; è necessario indicare il nome di ciascuno dei file di output. Le opzioni di Trimmomatic sono:

• PE indica che lavoriamo su file paired-end
• -trimlog indica che vogliamo salvare un file di log
• -threads indica il numero di processori che vogliamo utilizzare per lanciare il programma
• ILLUMINACLIP:ALL_ADAPTERS.fa:2:30:10 è la parte che indica la rimozione degli adattatori. Il file ALL_ADAPTERS.fa è il file in formato fasta che contiene la sequenza degli adattatori maggiormente utilizzati
• LEADING:16 indica la rimozione di N basi in testa alla sequenza, sotto la qualità indicata
• TRAILING:16 indica la rimozione di N basi in coda alla sequenza, sotto la qualità indicata
• SLIDINGWINDOW:4:25 effettua un trimming basato su una “finestra” di basi, rimuovendo le basi nel momento in cui la qualità media per quella finestra scende sotto la qualità indicata
• MINLEN:36 indica la rimozione delle reads che sono più corte della lunghezza indicata

# ci spostiamo nella cartella tc, creiamo nuova cartella e scarichiamo il programma
# usiamo percorso assoluto da qualsiasi punto del computer per spostarci dentro una qualsiasi cartella
cd /home/usr/Scaricati/tc # usr è specifico per il vostro utente (il mio è dlfptr70)
mkdir Trimmomatic
wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.36.zip -O Trimmomatic/Trimmomatic-0.36.zip # con l'opzione -O specifichiamo dove vogliamo scaricare il file e con che nome
# decomprimiamo senza cambiare cartella e rimuoviamo l'archivio
unzip Trimmomatic/Trimmomatic-0.36.zip -d Trimmomatic/ && rm Trimmomatic/Trimmomatic-0.36.zip
# controlliamo il funzionamento del programma
java -jar Trimmomatic/Trimmomatic-0.36/trimmomatic-0.36.jar -h # -h è l'opzione per help

A questo punto possiamo lanciare Trimmomatic sui nostri file; anche se dal controllo qualità non abbiamo identificato sequenze di adattatori come contaminanti, un classico approccio potrebbe essere quello di dare come input di file degli adattatori il file fornito dal programma che quindi rimuoverà tutte le sequenze corrispondenti a quelle che trova

# creiamo file unendo in uno tutti gli adattatori
cat Trimmomatic/Trimmomatic-0.36/adapters/*.fa > Trimmomatic/Trimmomatic-0.36/adapters/ALL_ADAPTERS.fa
# cat è un comando che unisce file "row-wise"; l'asterisco (*.fa) indica di selezionare tutti i file che hanno estensione .fa
java -jar Trimmomatic/Trimmomatic-0.36/trimmomatic-0.36.jar PE -threads 10 -trimlog \
Trimmomatic/log_trimmomatic \
Sample_Untreated_1_1_subset.fastq.gz Sample_Untreated_1_2_subset.fastq.gz \
Trimmomatic/Sample_Untreated_1_1_subset_trim_paired.fastq.gz Trimmomatic/Sample_Untreated_1_1_subset_trim_unpaired.fastq.gz \
Trimmomatic/Sample_Untreated_1_2_subset_trim_paired.fastq.gz Trimmomatic/Sample_Untreated_1_2_subset_trim_unpaired.fastq.gz \
ILLUMINACLIP:Trimmomatic/Trimmomatic-0.36/adapters/ALL_ADAPTERS.fa:2:30:10 \
LEADING:16 TRAILING:16 SLIDINGWINDOW:4:25 MINLEN:36 &> Trimmomatic/log_trimmomatic_screen
# con il simbolo &> andiamo a salvare su un file anche l'output che normalmente sarebbe stampato sullo schermo

Dovrebbe terminare dopo qualche secondo senza aver stampato a schermo alcuna informazione. Andiamo a vedere se ha creato i file ed esaminiamo i file di log (log_trimmomatic e log_trimmomatic_screen); il primo contiene le seguenti info per ciascuna read:

• the read name
• the surviving sequence length
• the location of the first surviving base, aka. the amount trimmed from the start
• the location of the last surviving base in the original read
• the amount trimmed from the end

Il secondo contiene le statistiche generali

# guardiamo prima le statistiche generali
less -S Trimmomatic/log_trimmomatic_screen
# poi l'altro file
less -S Trimmomatic/log_trimmomatic

A questo punto ripetiamo il controllo di qualità, questa volta sulle reads trimmate, solo sui file paired

# creiamo cartella dove salveremo nuovo controllo qualità
mkdir fastqc/results_after_trimming
# lanciamo FastQC sulle sequenze trimmate
./fastqc/FastQC/fastqc Trimmomatic/Sample_Untreated_1_1_subset_trim_paired.fastq.gz \
Trimmomatic/Sample_Untreated_1_2_subset_trim_paired.fastq.gz -o fastqc/results_after_trimming/

Andiamo ad esplorare i report e vedere se le statistiche sulla qualità sono cambiate

MAPPING

Ora i nostri file sono pronti per essere allineati sul genoma di riferimento (Arabidopsis). Per quasi tutti i programmi che effettuano allineamento il primo passaggio consiste della creazione dell’indice del genoma. Perché viene creato l’indice del genoma: citazione presa da qualcuno ma non ricordo chi (modificata) “Se entro in una biblioteca con la frase “La genetica è il ramo della biologia che si occupa del materiale ereditario” e voglio trovare il libro in cui è scritta, farò molta fatica a trovarlo, a meno che all’ingresso della biblioteca non si trovi un catalogo che mi indirizza almeno verso la sezione “Scienza””. L’indice quindi consiste in una serie di file che riassumono le principali caratteristiche del genoma di un organismo. Per effettuare l’allineamento useremo un programma chiamato STAR (https://github.com/alexdobin/STAR). Andiamo a scaricarlo insieme ai file aggiornati del genoma e del trascrittoma di Arabidopsis

# siamo nella cartella tc
mkdir STAR
# scarichiamo e salviamo nella cartella appena creata
wget https://github.com/alexdobin/STAR/archive/2.6.0c.tar.gz -O STAR/2.6.0c.tar.gz
# decomprimiamo e rimuoviamo archivio
tar -xf STAR/2.6.0c.tar.gz -C STAR/ && rm STAR/2.6.0c.tar.gz
# controlliamo il funzionamento
./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STAR -h

Per creare l’indice del genoma, abbiamo bisogno del genoma. Andiamo a scaricarlo dal sito di Ensembl Plants

# genoma
mkdir genome
wget ftp://ftp.ensemblgenomes.org/pub/release-39/plants/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz -O genome/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
# trascrittoma
wget ftp://ftp.ensemblgenomes.org/pub/release-39/plants/gtf/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.39.gtf.gz -O genome/Arabidopsis_thaliana.TAIR10.39.gtf.gz
# i file devono essere decompressi
gunzip genome/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
gunzip genome/Arabidopsis_thaliana.TAIR10.39.gtf.gz

A questo punto possiamo creare l’indice con il comando di STAR --runMode genomeGenerate. Vediamo nel dettaglio le altre opzioni:

• ./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STAR chiamiamo il programma
• --runThreadN numero di threads con cui vogliamo lanciare il programma
• --runMode indica che vogliamo creare l’indice
• --genomeDir indica la cartella dove vogliamo che venga salvato l’output
• --genomeFastaFiles cartella dove è presente il file del genoma
• --sjdbGTFfile indichiamo che vogliamo usare anche il trascrittoma per aiutarlo nella creazione dell’indice
• --sjdbOverhang specifica la lunghezza delle reads meno 1. Serve per creare il database delle giunzioni di splicing

Lanciamo il comando con le opzioni e i percorsi adattati al nostro caso

# creiamo cartella per l'indice all'interno della cartella del genoma
mkdir -p genome/index
# lanciamo la creazione dell'indice
./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STAR \
--runThreadN 10 \
--runMode genomeGenerate \
--genomeDir genome/index \
--genomeFastaFiles genome/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa \
--sjdbGTFfile genome/Arabidopsis_thaliana.TAIR10.39.gtf \
--sjdbOverhang 100

Con 10 core dovrebbe impiegare circa 3 minuti. Andiamo a vedere se ha creato i file nella cartella index

# spostiamo il file Log.out dove dovrebbe stare
mv Log.out genome/index/
# guardiamo i file appena creati
ls -lrht genome/index
# raccolgono diverse statistiche sul genoma di Arabidopsis, comprese info sul db di splicejunctions

Possiamo ora andare ad effettuare l’allineamento. Ovviamente per l’allineamento useremo le reads trimmate e pulite (le paired). Vediamo nel dettaglio i comandi:

• ./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STAR chiamiamo il programma
• --runThreadN numero di threads con cui vogliamo lanciare il programma
• --genomeDir dove sta la cartella con l’indice
• --readFilesIn qui scriviamo il percorso e il nome dei campioni da mappare
• --readFilesCommand qual è il formato di comrpessione dei file da mappare
• --outFileNamePrefix percorso e prefisso che vogliamo assegnare ai file di output
• --alignIntronMax lunghezza massima (stimata) degli introni del nostro organismo
• --outSAMtype tipo es estensione del file di output

Attenzione che la maggior parte dei programmi di allineamento sono creati con impostazioni di default basate sul genoma umano –> aggiustare i parametri sulla base del proprio organismo! Andiamo a lanciare il comando

# creiamo cartella dove salveremo file di allineamento
mkdir -p STAR/aligned
# lanciamo l'allineamento
./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STAR \
--runThreadN 10 \
--genomeDir genome/index/ \
--readFilesIn Trimmomatic/Sample_Untreated_1_1_subset_trim_paired.fastq.gz \
Trimmomatic/Sample_Untreated_1_2_subset_trim_paired.fastq.gz \
--readFilesCommand zcat \
--outFileNamePrefix STAR/aligned/SU_1_ \
--alignIntronMax 2000 \
--outSAMtype BAM SortedByCoordinate

Se è tutto a posto l’allineamento dovrebbe impiegare circa 10 secondi. Come vedete la lunghezza massima media per gli introni di Arabidopsis è di circa 2000 bp, e gli ho detto che voglio che il file di output sia in formato .bam (formato allineamento binario). Il file .bam è in un formato che non può essere aperto con il terminale, ma che andremo poi a visualizzare con un apposito software. Per completezza vi inserisco un immagine dello stesso file di allineamento ma in un formato visualizzabile e consultabile sul terminale. Il formato di questo file è .sam (Sequence Alignment Format), formato non molto utilizzato perché non essendo compresso occupa molto spazio

SAM file

La prima riga descrive il numero della versione del file

Dalla seconda alla ottava sono riportate le lunghezze in paia di basi dei cromosomi di Arabidopsis

La nona e la decima riga indicano il programma utilizzato per l’allineamento e il comando digitato

Dalla undicesima riga ogni riga coincide con una read, le colonne sono descritte di seguito

La colonna più interessante è forse la sesta (CIGAR) dove vengono specificati nel dettaglio gli allineamenti di ciascuna reads con il reference

SAM file description

Andiamo ora a vedere i nostri file creati

# il file che ci interessa è chiamato SU_1_Aligned.sortedByCoord.out.bam
ls -rht STAR/aligned/
# leggiamo qualche statistica in questo altro file
less -S STAR/aligned/SU_1_Log.final.out

Quanto vi risulta la percentuale di Uniquely mapped reads?

ALIGNMENT VISUALIZATION

Dopo aver effettuato l’allineamento può essere utile effettuare un controllo visivo sull reads allineate. A tal fine si utilizza spesso un software chiamato IGV (http://software.broadinstitute.org/software/igv/) che permette tra le altre cose di visualizzare i file dell’allineamento dopo aver caricato un genoma di riferimento

Integrative Genomics Viewer

Andiamo a scaricare il programma dal sito http://software.broadinstitute.org/software/igv/download cliccando in basso su “Launch with 2 GB”. Attendiamo lo scaricamento e apriamo il programma con doppio clic. Una volta aperto in alto a sinistra selezionare il genoma di A. thaliana TAIR10. A questo punto possiamo caricare i risultati del nostro allineamento, non prima però di aver creato un indice per il file .bam. Un altro indice direte voi, sì questa volta però non eseguiremo nessun passaggio ma vi farò prelevare il file dell’indice già preparato dal nostro server. Per completezza vi mostro il comando utilizzato per creare l’indice con Samtools

samtools index -b SU_1_Aligned.sortedByCoord.out.bam

Su questi computer non possiamo installare il programma che si utilizza, tra le altre cose, anche per la creazione degli indici per i file di allineamento. Il programma si chiama Samtools (http://www.htslib.org/) ed è molto utilizzato in bioinformatica perché svolge numerose funzioni. Andiamo a copiare nella cartella dove sta il file .bam appena creato da noi, il file dell’indice, che permette che il file .bam venga visualizzato correttamente in IGV

scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/SU_1_Aligned.sortedByCoord.out.bam.bai STAR/aligned/
# il file che vedete deve avere lo stesso nome del file di allineamento, ma ha un'estensione diversa, .bai

A questo punto potete tornare su IGV e cliccare su “File -> Load from file” e selezionare dalla cartella aligned il file SU_1_Aligned.sortedByCoord.out.bam. Selezioniamo Chr1 e zoomiamo sino a che la rotella inizia a girare, a quel punto possiamo visualizzare le nostre reads allineate ai geni del genoma di riferimento

Non sarà così facile trovare un po’ di reads!

Nel caso non funzionasse, copiare dal server alla vostra cartella, sia il file .bam che il .bai (potrebbe essere che il bai debba essere creato dallo stesso bam per funzionare la visualizzazione)

# bam
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/SU_1_Aligned.sortedByCoord.out.bam STAR/aligned/
# bai
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/SU_1_Aligned.sortedByCoord.out.bam.bai STAR/aligned/

La lunghezza degli introni è stata rispettata?

READS COUNTING

Una volta controllato che l’allineamento è andato a buon fine e siamo soddisfatti di come il programma di mapping ha lavorato, possiamo spostarci nel prossimo step verso l’analisi dei geni, ovvero la conta delle reads. Per conta delle reads si intende il contare quante reads sono state assegnate (hanno “mappato”) a un determinato gene (o trascritto, o esone, le cosiddette genomic features) in seguito al mapping. Per questo passaggio utilizzeremo il software featureCounts (http://subread.sourceforge.net/) del pacchetto Subread. Scarichiamo e testiamo il programma

# nella cartella tc/ creiamo cartella
mkdir featureCounts
# scarichiamo dal server perché più comodo
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/featureCounts featureCounts/
# testiamo funzionamento
./featureCounts/featureCounts -v

Vediamo le opzioni e i parametri di featureCounts:

• featureCounts lanciamo il programma
• -a forniamo il file di annotazione del trascrittoma
• -o percorso e nome file di output
• -p indichiamo che lavoriamo su reads paired e che verranno contati i frammenti e non le reads
• -t indica quale feature contare
• -g indica che vogliamo raggrppare tutte le feature in questa meta-feature
• -T quanti threads vogliamo usare per correre il programma
• in.bam nome del file di allineamento che vogliamo come input

A questo punto andiamo a lanciare il programma, ricordandoci che nel nostro caso vogliamo contare le feature del file di allineamento .bam

# creiamo cartella per output
mkdir -p featureCounts/output
# lanciamo la conta
./featureCounts/featureCounts \
-a genome/Arabidopsis_thaliana.TAIR10.39.gtf \
-o featureCounts/output/SU_1_counts.txt \
-p \
-t exon \
-g gene_id \
-T 10 \
STAR/aligned/SU_1_Aligned.sortedByCoord.out.bam

Dovrebbe terminare in qualche secondo. Andiamo ad esplorare i file di output

ls -lrht featureCounts/output/
less -S featureCounts/output/SU_1_counts.txt.summary
less -S featureCounts/output/SU_1_counts.txt
# vediamo solo le info che ci interessano
cut -f 1,7 featureCounts/output/SU_1_counts.txt | sed 1d | less -S

Le colonne che ci interessano sono la prima e l’ultima (7). Perché per alcuni geni abbiamo più di un valore per le colonne 2, 3, 4 e 5? I passaggi che abbiamo visto finora in una situazione normale andrebbero ripetuti per tutti i campioni che abbiamo fatto sequenziare, sino ad ottenere una matrice di conte, come l’ultimo file che abbiamo visto, per ciascun campione. Da questa matrice di conte (counts matrix) creeremo un file unico che sarà l’input per la fase successiva. Visto che nella nostra esercitazione abbiamo elaborato solo un file, vi fornirò io la counts matrix completa (6 campioni, 3 controlli e 3 trattati)

# scarichiamo dal server il file con le conte per i 3 controlli e i 3 trattati
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/countDF_unvs3h featureCounts/
# vediamo il file
less -S featureCounts/countDF_unvs3h

INSTALLING R AND RSTUDIO

A questo punto abbiamo concluso la prima fase del corso che abbiamo svolto usando solamente il terminale. Per la seconda fase utilizzeremo invece il software R e il suo ambiente di sviluppo integrato (IDE) chiamato RStudio. Per prima cosa vi invito a scaricare e installare R ed RStudio seguendo le istruzioni di seguito a seconda del vostro SO:

• Windows: cliccare su http://cran.mirror.garr.it/mirrors/CRAN/bin/windows/base/R-3.5.0-win.exe, attendere lo scaricamento e installare il file con estensione .exe appena scaricato, lasciare le impostazioni di default. Per installare RStudio andare alla pagina https://www.rstudio.com/products/rstudio/download/#download e alla voce “Installers for Supported Platforms” cliccare su “RStudio 1.1.453 - Windows Vista/7/8/10”, attendere lo scaricamento e installare il file con estensione .exe appena scaricato, lasciare le impostazioni di default.

• MacOS: cliccare su http://cran.mirror.garr.it/mirrors/CRAN/bin/macosx/R-3.5.0.pkg, attendere lo scaricamento e installare il file con estensione .pkg appena scaricato, lasciare le impostazioni di default. Per installare RStudio andare alla pagina https://www.rstudio.com/products/rstudio/download/#download e alla voce “Installers for Supported Platforms” cliccare su “RStudio 1.1.453 - Mac OS X 10.6+ (64-bit)”, attendere lo scaricamento e installare il file con estensione .dmg appena scaricato, lasciare le impostazioni di default.

• Linux: informazioni disponibili all’indirizzo https://medium.com/@GalarnykMichael/install-r-and-rstudio-on-ubuntu-12-04-14-04-16-04-b6b3107f7779

Per qualsiasi problema scrivetemi all’indirizzo delfino.pietro@gmail.com o possiamo vedere domani l’installazione in aula, ma sarebbe meglio arrivare già con entrambi i programmi installati e alcuni pacchetti che elenco di seguito, con i rispettivi comandi per l’installazione. Aprire RStudio e nella console digitare i seguenti comandi con sfondo grigio, in ordine, uno per uno:

• Bioconductor
  • source("https://bioconductor.org/biocLite.R")
  • biocLite()
• tidyverse
  • install.packages("tidyverse")
• DESeq2
  • source("https://bioconductor.org/biocLite.R")
  • biocLite("DESeq2")
• edgeR
  • source("https://bioconductor.org/biocLite.R")
  • biocLite("edgeR")
• biomaRt
  • source("https://bioconductor.org/biocLite.R")
  • biocLite("biomaRt")
• clusterProfiler
  • source("https://bioconductor.org/biocLite.R")
  • biocLite("clusterProfiler")

DIFFERENTIAL GENE EXPRESSION ANALYSIS

setwd("/home/usr/Scaricati/tc/featureCounts") 
# "usr" è variabile a seconda del vostro nome utente, usate il tab per completare il percorso

# carichiamo la libreria
library(DESeq2)

# carichiamo il file
cts <- as.matrix(read.table("featureCounts/countDF_unvs3h", sep="\t", header=T, col.names = c("ID","un1","un2","un3","NO3_1","NO3_2","NO3_3"), row.names = "ID"))

# non eseguire questo comando, è solo per me
cts <- as.matrix(read.table("Scaricati/tc/featureCounts/countDF_unvs3h", sep="\t",
                            header=T, col.names = c("ID","un1","un2","un3","NO3_1","NO3_2","NO3_3"),
                            row.names = "ID"))
# vediamo cosa abbiamo caricato
head(cts)

##           un1  un2  un3 NO3_1 NO3_2 NO3_3
## AT1G01010 248  535  190   830   839   497
## AT1G01020 255  684  943   842   679   909
## AT1G01030 546  607  568    61    45    84
## AT1G01040 690 1970 1386   607   687   419
## AT1G01046   2    7    7     1     3     0
## AT1G01050 597 1252 2008  1015  1054  1061

# creiamo oggetto per il disegno sperimentale
coldata <- data.frame(row.names = colnames(cts))
coldata$condition <- c(rep("untreated", 3), rep("treated", 3))
coldata$time <- c(rep("0", 3), rep("3", 3))
str(coldata)

## 'data.frame':    6 obs. of  2 variables:
##  $ condition: chr  "untreated" "untreated" "untreated" "treated" ...
##  $ time     : chr  "0" "0" "0" "3" ...

# convertiamo le colonna in factor per poterle usare nella formula
coldata$condition <- as.factor(coldata$condition)
coldata$time <- as.factor(coldata$time)
str(coldata)

## 'data.frame':    6 obs. of  2 variables:
##  $ condition: Factor w/ 2 levels "treated","untreated": 2 2 2 1 1 1
##  $ time     : Factor w/ 2 levels "0","3": 1 1 1 2 2 2

# costruiamo l'oggetto DESeq
dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData = coldata,
                              design = ~ condition )

# ordiniamo i livelli per usare "untreated" come confronto
dds$condition <- relevel(dds$condition, ref = "untreated")

# vediamo info sull'oggetto creato
dds

## class: DESeqDataSet 
## dim: 33602 6 
## metadata(1): version
## assays(1): counts
## rownames(33602): AT1G01010 AT1G01020 ... ATMG01400 ATMG01410
## rowData names(0):
## colnames(6): un1 un2 ... NO3_2 NO3_3
## colData names(2): condition time

# rimuoviamo geni mai espressi
dds <- dds [rowSums(counts(dds)) > 0, ]

dds

## class: DESeqDataSet 
## dim: 24383 6 
## metadata(1): version
## assays(1): counts
## rownames(24383): AT1G01010 AT1G01020 ... ATMG01400 ATMG01410
## rowData names(0):
## colnames(6): un1 un2 ... NO3_2 NO3_3
## colData names(2): condition time

# quanti erano?